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Abstract. Point-source contamination in high-precision Cosmic Microwave Background (CMB) maps severely 
affects the precision of cosmological parameter estimates. Among the methods tha t have been prop osed for source 
detection, wavelet techniques based on "optimal" filters have been proposed (e.g., sanz et al. 2001). In this paper 



we show that these filters are in fact only restrictive cases of a more general class of matched filters that optimize 
signal-to-noise ratio and that have, in general, better source detection capabilities, especially for lower amplitude 
sources. These conclusions are confirmed by some numerical experiments. 
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1. Introduction 

The separation of different physical components is an 
important issue in the analysis of Cosmic Microwave 
Background (CMB) data. Among the foreground compo- 
nents, point-sources deserve especial attention given their 
extreme non-Gaussianity and highly variable spectral in- 
dex. A brief summary of different methods can be found 



in 3anz et al. (2001) (henceforth SHAf) 



Methods used to detect point sources should act as 
high-pass filters to detect the high frequency structure in- 
troduced by the sources in the CMB data while masking 
other foreground contamination (like dust, synchrotron 
and free-free emission) characterized by lower frequen- 
cies. In principle, the good space/frequency characteris- 
tics of wavelets should make these functions an attractive 



nothing seems to be gained with these wavelet filters and 
that other simpler techniques lead to better source detec- 
tion methods. 

We first set up the basic framework. Although the 
detection of point-like sources in CMB maps is a two- 
dimensional problem, we present our arguments in R", as 
in SHM, because the same methods may be used in other 
applications. 

The sources are assumed to be point-like signals con- 
volved with the beam of the measuring instrument and 
are thus assumed to have a known profile t{x). The sig- 
nal y{x), X S i?", is modeled as 



■i{x) 



choice for such tasks (e.g., Cayon et al. 2000; Vielva et al 
2001a| ,[b|, and references therein). Indeed, wavelets have 
been proved optimal for detection of point like singulari- 
ties - at least for one-dimensional signals. Point sources in 
CMB maps, however, are not point singularities because 
the signal is smoothed by the beam of the antenna; sources 
are expected to have the shape of the antenna's pattern. 
The question is then how to include this beam profile in- 
formation in the wavelet analysis. SHM considered this 
question and, were lead to define optimal pseudo-filters 
for source detection in CMB maps. However, we show that 



where 



SjN ^Ajt{x~Xj), 



(1) 



(2) 



Aj and Xj are, respectively, unknown source amplitudes 
and locations, and z{x) is a zero-mean background with 
power-spectrum P{q) 



E[ziq)z*iq')] = P{q) 5"(g-g')- 



(3) 
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Henceforth E[-] and " * " will denote the expectation and 
complex conjugate operators, respectively, S^{q — q') the 
n-dimensional Dirac distribution, and z(q) the Fourier 
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transform of z{x) 



4q) = f 



z{x) e 



dx. 



(4) 



To properly remove the point sources from the signal 
we need to estimate the locations and amplitudes 

(fluxes) {Aj\ of the sources. 

A classical method used to estimate source locations 
is based on identifying peaks in the cross-correlation func- 
tion 

/ + 00 
y{x + h) r(6) dh. (5) 
-OO 

The rationale is that c(a;o) measures the similarity be- 
tween the source profile with a section of y(a^) centered at 
peak in c{xo) is an indication of a source signal at 
Xo- Eq. (|^) is a filtering of the signal yix) with a filter r(b) 
that amplifies the characteristic frequencies of the source. 
Once the sources have been located, their amplitudes can 
be estimated by means of classical fitting procedures like 
least squares. 

The cross-correlation function (H) docs not take into 
account the background characteristics. This is a great 
disadvantage in cases where the power spectrum P{q) is 
known or a good estimate is available. In Sects, y and 
^ we consider other methods that take into account this 
information and that may be considered extensions of the 
cross-correlation method. 

The basic procedure we consider is as follows. The sig- 
nal is first filtered to enhance the sources with respect 
to the background. This is done by cross-correlating the 
signal y(x) with a filter $ as in (with $ in place of 
r). The source locations are then determined by selecting 
the peaks in the filtered signal that are above a selected 
threshold. Finally, the source amplitudes are estimated 
with the values of the filtered signal at the estimated lo- 
cations. The question we consider first is the selection of 
an optimal filter $ for such procedure. 

2. Designing an optimal filter 

The optimality criteria we use are based on the following 
assumptions (for futher justification of these assumptions 
see SHM). The source profile and background spectrum 
are known. The profile is spherically symmetric, charac- 
terized by a scale i?s, and the background is assumed to 
be isotropic. In this case, we can write s{x) = s{x), where 
X = \\x\\, and P{q) = P{q) for q = \\q\\. In addition, source 
overlap is assumed negligible. 

We consider the general family of spherically symmet- 
ric filters ^{x;b) of the form 

^x;b)=^{\\x-b\\) (6) 

with Fourier transform 4>{q). The filtered field is 



y{x) ^{x\ b) dx 



+ 00 



y{q) m e*-?" dq, 



(7) 



where y(<7) and (/)(g) are, respectively, the Fourier trans- 
forms of y{x) and (^{x). The mean and variance of w{b; (j)) 
are 



H{b; 



Ew{b](t>) 

r + co 

a / 

Jo 



s{q)4>{q)^"''''dq 



a^i<t>)^J^w^ib,<^>)-^,{b;<t>f 

r+oo 

= a/ q'^-' P{q)<Piq)dq, 
Jo 



(8) 



(9) 



where a = r-i(n/2). 

The first constraint on the filter concerns the second 
stage of the procedure; the source locations are assumed 
known and the objective is to estimate the amplitudes. 
Given the assumed distance between sources, it is enough 
to consider a field y{x) as in (|^) with a single source at the 
origin, s{x) — At{x). To estimate its amplitude we ask 
that w{0; (p) be an unbiased estimator of A - i.e., /i(0; cj)) — 
A - so that (j) is required to satisfy the equation 



7"-V(g)0(q) dq 



1 



(10) 



To enhance the magnitude of the source relative to the 
background we determine the filter $ that minimizes the 
variance cr^((/)). This has the effect of maximizing, among 
unbiased estimators, the detection level 



A^(O;0: 

a(</>) 



(11) 



which measures the capability of the filter to correctly 
detect a source at the prescribed location (see SHM). 

Since is chosen so that w{Q] cj)) is a minimum variance 
linear - in y{x) - unbiased estimator of A, it follows that 
(Gauss-Markov theorem) w{0; (p) is the (generalized) least 
squares estimate of A achieved by the filter 



1 riq) 



aa P{q) ' 
with minimum variance 



a^(0) 



1 

aa 



dq, (12) 



(13) 



Filter (|lj) is a particular case of a well known class of 
filters, known as matched filters in the engineering liter- 
ature, that are designed to optimize signal-to-noise ratio 
(e.g., iKozma fc Kelley 1965| ; |Pratt 199l] ). The argu- 



ments in this section also show that the filter introduced 



by Tegmark fc Oliveira-Gosta (1998 ) is a particular case 
of the matched filter (|l2|) with P{q) representing the back- 
ground power spectrum before the smoothing of the an- 
tenna. 

For white noise, P(q) — const — D, filter ( p^ simpli- 
fies to 

0(.) = ^ riq) (14) 
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which, up to a constant factor, is identical to the filter 
used in the classical cross-correlation function (^. This 
provides a justification for the use of the cross-correlation 
for source detection in white noise. 

Having recognized w{0;4>) as a least squares estima- 
tor of A, we close this section with some remarks from 
least squares methodology that we consider relevant. First 
note that, regardless of the spectrum P{q), the source pro- 
file T{q) properly normalized, that is (^(g) — T{q)/K for 
K = a J q'^~^T{q) (j){q) e^'^ ^'dg, also provides an unbiased 
estimator w(0; </>) of A. This is the (ordinary) least squares 
estimate that does not take into account the covariance 
of the background; it is unbiased but not minimum vari- 
ance. However, it is well known that when the covariance 
is actually estimated from the data, the ordinary least 
squares estimate may be better than the generalized one 
(e.g.. Draper fc Smith 1998 ). In other words, uncertaini- 



where A = ac — 6^ , 



ties in the spectrum estimates may lead to worse ampli- 
tude estimates than those obtained with the simpler cross- 
correlation filter. Uncertainties in the spectrum will also 
affect the selection of a detection threshold. 

Note also that the unbiasesness of w{0; </>) as an estima- 
tor of A depends on knowing the correct source location, 
it is not necessarily unbiased once the source locations are 
estimated from the data. This is shown in Sect. 



4.2 



3. Pseudo-filters 

In the pseudo-filter approach of SHM the filters are of the 
form with an additional scale dependence 



^{x;R,b) = 



\x - b\ 
R 



(15) 



for some spherically symmetric function ip. From (|15|) 
we get a set of wavelets of different scales R and shifts 
b provided the regularity conditions T{q)/P{q) 0, 
(l/P){dT/dhiq) ^ for 9 ^ are satisfied. The filtered 
field w{b, R; iji) at scale R is defined as in (^ but with 
ip{Rq) in place of 4>{q)- 

To determine an optimal filter ^, SHM minimize the 
variance of the filtered field subject to two constraints: 
First, w(0, i?o; ip) is required to be, as in the previous sec- 
tion, an unbiased estimator of A for some known i?o ~ -Rs ■ 
For the second constraint ip is selected so that ^(0, R; tp) 
has a local maximum at scale i?o- This constraint trans- 
lates to 



+ 00 



(7"-V(g)^(i?og) (n 



dlnr 
ding 



dq = 0. (16) 



Minimizing a'^(RQ;ip) with the two constraints yields the 
fihcr {SHM) 



i^iRoq) = 



1 r(g) 



a A P{q) 



nb + c — {na + b) 



dhiT{q) 
dhiq 



(17) 



T dr 



'"io Pdlnq^ 



(18) 



dq, 



la P \d\nq ^ 

and a is as in (^2|). This filter provides a field of variance 

n^a + 2nb + c 



(7^(i?o;V) 



tA 



(19) 



The estimator of the amplitude A obtained with this 
filter is again linear and unbiased and therefore, by the op- 
timality of least squares, a'^{Ro;ip) > cr^(0) regardless of 
the source profile and background spectrum. This means 
that the detection level of <& is at least as high, or higher, 
than that achieved with '5. 

3.1. /s ^ optimal for source detection? 

We have seen that the second constraint (|6|) does not 
increase the detection level when the source locations are 
known. But this is not surprising since the constraint is 
defined to take advantage of the known source scale to 
help determine source locations. We will show that even 
when source location uncertainty is taken into account, 
(|l6|) does not improve on the simpler filter $ based on the 
single contraint (px|). In other words, enough information 
about the scale of the source is already included in the 
derivation of the matched filter. Therefore, it seems that 
nothing is gained with wavelet methods. 

In principle, as explained by SHM, constraint (|6|) can 
be used to test if a detection corresponds to a true source 
by checking for its maximum at scale Rq. But finding a 
spike at the correct scale is not enough to make sure it is 
not a spurious noise artifact, we have to follow it across 
scales to make sure that it scales appropriately. In any 
similar test can also be designed for the matched 
filter: given that the source profile is assumed known, the 
matched filter leads to an amplitude versus scale depen- 
dence (see example in Sect. 4T) that can be determined 
and used for such a test. This functional dependence pro- 
vides more information than the existence of a local max- 
ima and should, in the presence of noise, lead to better 
location estimates. 



4. Examples 

To compare the theoretical performances of the filters $ 
and ^' we use the gain as defined by SHM 



9i'4',(p) 



V{^)/V{ij,Ro) 
a(i?o;^)/cT(0), 



(20) 



where 'D{iP,Rq) is as in (|ll|) but defined for ^{x;Ro,0). 
We start with a specific example with Gaussian profiles 
and power-law spectra. 
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Y=0.5 




Fig. 1. Filters '^{x; Rq,0) and ^{x) for n = 1 and white 
noise. Panel (a) shows the filters in spatial domain and 
panel (b) shows their power- spectra. 





7 = 2 





10 15 



Fig. 2. Cumulative energies and i?^ of ^'(a;; Rq, 0) and 
^{x), respectively, for n — 1 and different values of 7 (see 
text). 



4.1. Gaussian sources and P{q) = Dq 
Gaussian sources 



T(g) =re-(''«)'/^ 



(21) 



where 9 is the "standard deviation" defining the scale, 
provide an important family of source profiles. Indeed, 
in many practical applications the instrument's profile, 
and thus the point sources, are effectively characterized 
by Gaussian profiles. For the noise process we take the 
power-law spectrum P{q) = Dq~'^ . This family of spectra 
can be used to locally approximate spectra of other homo- 
geneous processes and has as special cases white (7 = 0) 
and 1// (7 = 1) noise processes. 

For a Gaussian profile and a power-law spectrum, Eq. 
(|l|) leads to 



_ r(n/2) {qoy 



r(m) 



r"/2 



(22) 





Fig. 3. Same as Fig. g but with n — 2. 

where m = {n + 7)/2. ^From (|3|) and ([l9| ) we obtain the 
following variance and gain 



D 



r(»/2) 

1/2 



1 



{n-1? 
4m 



(23) 



We see that g > 1, which shows that (/) has a higher 
detection level than ijj. For n — j — for example, one- 
dimensional process with 1// noise — the two filters lead 
to the same detection levels. This is expected since (l){q) = 
ip{Rq) for n = 7 (compare Eq. (|2|) with Eq. (24) in SHM). 
In other words, the second constraint is just redundant in 
this case. Note also that (j){q) is the Mexican hat wavelet 
for n = 2 and 7 = 2. This results sh ow that, contrary t o 
what has been claimed before (e.g., Cayon et al. 2000 ), 
the optimality of the Mexican hat wavelet does depend on 
the background spectrum. 

4.2. A numerical experiment 

We have presented theoretical arguments showing that 
(/)(g) has a better source detection capability than ^j{Rq) 
when the source location is known. We now confirm that cf) 
is still better when the uncertainity of source location es- 
timates is taken into account. To compare with the results 
in SHM, it is enough to consider a simple example with 
one-dimensional Gaussian sources and white noise (7 = 0) 
(we already know that (jj = for n = 1 and 1// noise). In 
this case (|l^ ) and (|3|) become, respectively. 



(t>{q) = ^ e 



-{lSf/2 



1/2 



(24) 



> 1. 



That is, the detection level of (p is about 20% larger than 
that of tp. The amplitude dependence on of a Gaussian 
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Filter 




S/N = 


1 




S/N = 


2 




S/N = 


3 




iV, 


A 


Nc 




A 






A 




23 


10 


1.83 


87 


10 


1.10 


100 


10 


0.97 




(15) 


(5) 


(1.94) 


(80) 


(5) 


(1.14) 


(100) 


(5) 


(0.98) 




31 


5 


1.55 


97 


5 


1.02 


100 


5 


0.98 



Table 1. Results of the numerical simulations concerning the detection capabilities of the matched and the pseudo- 
filters (see text) in the one dimensional case and for a white noise background. Nc, Ni are, respectively, the average 
number of correct and incorrect detections. Five hundred simulations were carried out using 100 gaussian sources of 
amplitude A — 1 and scale 9 — 3, regularly distributed along an array of about 16000 elements. The 3a detection level 
has been determined by filtering, via <i> and independent realizations of the background process. The amplitude 
estimates are based on the average amplitudes of all detected sources. The S/N is defined by A/ae where Ce is the 
standard deviation of the noise process. The second row for corresponds to results with a threshold (> 3a) chosen 
to achieve the same average number of incorrect detections as $. 





S/N 


= 1 


S/N 


= 2 


S/N 


= 3 


Filter 




A 


Pc 


A 


Pc 


A 




0.89 


1.14 


1.00 


1.01 


1.00 


1.00 




0.99 


1.03 


1.00 


1.00 


1.00 


1.00 



Table 2. Results of two-dimensional simulations with a 128 x 128 grid and white noise. A is as in Table I. Pc is the 
estimated probability of correctly detecting a source keeping the rate of incorrect detection at the same level as in 
Table ui Here ^' coincides with the Mexican hat wavelet. 



Filter 


S/N 


= 1 


S/N 


= 2 


S/N 


= 3 


Pc 


A 


Pc 


A 


Pc 


A 


MH 


0.38 


1.79 


0.93 


1.11 


1.00 


1.03 




0.40 


1.78 


0.93 


1.11 


1.00 


1.03 




0.47 


1.60 


0.97 


1.07 


1.00 


1.02 



Table 3. Results of two-dimensional simulations with a 128 x 128 grid and 1// noise. A is as in Table |. Pc is the 
estimated probability of correctly detecting a source keeping the rate of incorrect detection at the same level as in 
Table |l| MH stands for Mexican hat wavelet. 

filters are quite different. For example, to provide filtered 
sources with the scale Rq , has to pass higher frequencies 
than (j). This can be a problem for signals contaminated 
by high frequency noise. 

Table ^ shows the average number of correct and in- 
correct detections obtained with and <& and a fixed 3a 
threshold for signal-to- noise (S/N) ratios equal to 1, 2 and 
3. The two filters give equivalent results for higher S/N 
sources. We see that (j) leads to a higher number of cor- 
rect detections and a lower number of incorrect ones. But 



source of scale Rs filtered with (f> is 

A fit to this dependence can be used to determine if a de- 
tection corresponds to a source of the correct scale, just as 
large wavelet coefficients would be tracked across different 
wavelet scales. 

Fig. |l| shows the filters ^{x) and '^{x; Ro, 0) and their 
corresponding Fourier transforms. It shows that the two 



6 



R. Vio, L. Tenorio & W. Wamsteker: Optimal Detection of Sources 



a proper comparison should take into account that the 
filters require different thresholds. The second row for 4" 
shows the corresponding results when the threshold is cho- 
sen to lead to the same average number of rejections as $. 
For low S/N sources $ leads again to a higher number of 
correct detections while for larger S/N they give similar 
results. 

To compare amplitude estimates that include location 
uncertainty, we take the average of the amplitudes (since 
all the generated sources have the same amplitude) of all 
detections. The results are shown in Table |l|. The errors 
in the amplitudes are of the order of 0.5% or less. We see 
that amplitude estimates are biased when the source lo- 
cations are estimated, and that the bias is larger for ^P. 
For low S/N the amplitude is overestimated because high 
peaks are easier to detect and noise peaks are incorrectly 
classified as sources. For high S/N sources we also have 
centering problems but this time the smaller amplitude 
in the noise peaks leads to underestimated source ampli- 
tudes. We can draw similar conclusions from the results 
of two-dimensional simulations shown in Tables ^ and |^; 

To conclude, note that the spatial support of the filters 
is an important factor when the assumption of nonover- 
lapping filters is invalid. The support must be small com- 
pared to the distance between the sources. Figs. ||-^ show 
(for 71 = 1, 2) the cumulative energies 



and 



E^(.t) 



E$(a;) 



/;^^(6;i?o,0)d& 
j;;°^>^b;Ro,0)db 



(26) 



(27) 



as a function of 7 for a background spectrum P{q) = 
Dq~^ . These functions measure the "energy concentra- 
tion" of the filters and provide information about their 
"spatial" support. The figures show that also in this re- 
spect $(x) has similar, if not better, characteristics than 
^'(a;; Rq, 0). In particular, the filter $ has a tighter support 
for faster decaying noise spectra. 

4.2.1. Selecting the detection level 

We briefly justify our selection of 3ct level. This threshold 
should be chosen large enough to reduce the number of 
false detections but small enough not to miss too many 
sources. To properly choose a threshold we have to un- 
derstand the statistics of local maxima above a chosen 
level. For a general homogeneous (Gaussian) random field 
this is a difficult question (some asymptotic results can be 
found in Adler 1981) but simulations can be carried out 
when the background power spectrum is known. Our sim- 
ulations showed that the traditional 5cr level is too con- 
servative for the signal lengths used in the examples. If 
Pm is the proportion of local maxima above 3(t for a field 
without sources. We found that the probability that Pm is 
higher than 0.001 is about 9% for the signal filtered with 
\I' and less than 10"^ for the signal filtered with 



5. Summary and Conclusions 

We have revisited the problem of estimating point sources 
of known profile in an isotropic background. The methods 
we considered are based on two basic interrelated proce- 
dures: source detection by thresholding of local maxima, 
and amplitude estimation by linear filtering. We have com- 
pared the effects of different constraints on the selection 
of an optimal filter. 

The first constraint is typical in matched filter method- 
ology where S/N is maximized. The optimal filter pro- 
vides unbiased least squares estimates of source ampli- 
tudes at known source locations. By the optimality of least 
squares, these amplitude estimates can not be improved 
with any other unbiased linear filter. However, uncertaini- 
ties in source locations introduce a bias in amplitude esti- 
mates. Amplitudes are overestimated at low S/N and un- 
derestimated at high S/N. A second constraint introduced 
by SHM is designed to improve estimates of source loca- 
tions by maximizing amplitudes at the correct scale. We 
found that this constraint does not lead to better estimates 
as compared to those obtained with a simpler matched fil- 
ter, especially for lowe r S/N sources. For hi gh S/N sources, 
as it was observed by Cayon et al. (2000 ), the results of 
the two methods are the same. These results contradict 
previous optimality studies of wavelet-based filters for de- 
tection of point-sources of known profile in an isotropic 
CMB background of known spectrum. 
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